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ABSTRACT 

We apply our two-dimensional (2D), radially self-similar steady-state aeeretion 
flow model to the analysis of hydrodynamic simulation results of supercritical accre¬ 
tion flows. Self-similarity is checked and the input parameters for the model calcula¬ 
tion, such as advective factor and heat capacity ratio, are obtained from time-averaged 
simulation data. Solutions of the model are then calculated and compared with the 
simulation results. We find that in the converged region of the simulation, excluding 
the part too close to the black hole, the radial distribution of azimuthal velocity v^, 
density p and pressure p basically follows the self-similar assumptions, i.e. they are 
roughly proportional to r“", and respectively, where n ~ 0.85 for the mass 

injection rate of 1000Le/c^, and n ~ 0.74 for 3000Le/c^. The distribution of Vr and 
V 0 agrees less with self-similarity, possibly due to convective motions in the rQ plane. 

The distribution of velocity, density and pressure in Q direction obtained by the steady 
model agrees well with the simulation results within the calculation boundary of the 
steady model. Outward mass flux in the simulations is overall directed toward polar 
angle of 0.8382 rad (~ 48.0°) for 1000Le/c^, and 0.7852 rad (~ 43.4°) for 3000Le/c^, 
and ~94% of the mass inflow are driven away as outflow, while outward momentum 
and energy fluxes are focused around the polar axis. Part of these fluxes lie in the re¬ 
gion that are not calculated by the steady model, and special attention should be paid 
when the model is applied. 
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1. Introduction 


Recent development in observations has shown growing evidence that a large portion of mass 
is blown away from accretion flow onto black holes in the form of outflows. For example, the 
accretion rate onto the Galactic Center is estimated to be about 10'^ Mq yr“^ (Baganoff et al. 
2003) at the outer boundary, while the detected high linear polarization at radio waveband limits the 
mass inflow rate near the event horizon to be < 10“^ Mq yr“^ (Baganoff et al. 2003) or 10“’ - 10“^ 
Mq yr“^ (Marrone et al. 2007), which implies that most of the accreting material cannot reach the 
black hole. Outflows are also observed through blue-shifted absorption lines of galactic sources 
(e.g.. Miller et al. 2004; Kotani et al. 2006; Kubota et al. 2007; Neilsen et al. 2012; Ponti et al. 
2012) and active galactic nuclei (AGNs; Terashima & Wilson 2001; Pounds et al. 2003; Reeves et 
al. 2003; Ganguly & Brotherton 2008; Pounds & Reeves 2009). 

Theoretically, it is believed that outflows are likely to be generated when the mass supply rate 
is much less than (’underfed’) or much greater than (’overfed’) the Eddington accretion rate. The 
’underfed’ case can be described by an optically thin advection-dominated accretion flow (ADAF; 
Narayan & Yi 1994, 1995a, b; Abramowicz et al. 1995; for reviews, see Narayan et al. 1998; Kato 
et al 2008; Narayan & McClintock 2008). In an ADAF, energy released via viscous dissipation 
cannot be radiated away efficiently, which causes the gas to be heated to very high temperature 
and become unbound, and strong outflows are likely to be generated (e.g. Narayan & Yi 1994; 
Quataert & Narayan 1999). The ’overfed’ case can be described by a ’slim disk’ (Abramowicz 
et al. 1988; Beloborodov 1998; Chen & Wang 2004; S^dowski 2009; Takeuchi et al. 2009; 
etc.). While the original slim disk model does not account for outflows due to treatment limit, it 
has been discovered that there is an upper limit of mass supply rate for slim disks beyond which 
outflows become inevitable (Gu & Fu 2007; Jiao et al. 2009; Jiao & Fu 2009; Gu 2012, 2015). 
In fact, even for accretion rates below this limit, outflows are likely to be generated due to strong 
radiation pressure (e.g. Ohsuga et al. 2009). In both cases, the existence of outflows has been 
validated by numerical simulations (e.g.. Stone et al. 1999; Igumenshchev & Abramowicz 2000; 
Stone & Pringle 2001; Hawley et al. 2001; Hawley & Balbus 2002; McKinney & Gammie 2002; 
Igumenshchev et al. 2003; De Villiers et al. 2003; Okuda et al. 2005; Ohsuga et al. 2005, 2009; 
Ohsuga & Mineshige 2007, 2011; Yuan & Bu 2010; McKinney et al. 2012; Narayan et al. 2012; 
Yuan et al. 2012a, 2012b; Takeuchi et al. 2013; etc.). 

However, it is very difficult to apply numerical simulation results directly to observations, 
e.g. fitting the spectra of accretion-powered astrophysical systems, due to the heavy computation 
they require. Simple analytic steady models are still the only accessible way of making direct 
link between theory and observations for astronomers. Nevertheless, steady models which contain 
outflows are quite limited, among which are contributions from Narayan & Yi (1995a), Xu & Chen 
(1997), Blandford & Begelman (1999, 2004; hereafter BB99, BB04), Xue & Wang (2005), Gu et 



-3 


al.(2009), S^dowski et al.(2011), Jiao & Wu (2011; hereafter JWll), Begelman (2012), etc. (see 
JW11 for a detailed discussion of the contributions and caveats of some of these works). Although 
these models are established self-consistently, it is not clear whether the assumptions and results 
of them agree with the simulation results, which in some sense represent the structure of real 
astrophysical accretion flows. 

In this paper, we contribute to this topic by applying our 2D self-similar accretion flow model 
(JWll) to the analysis of simulation results. As the first step, we select two samples of 2D 
radiation-hydrodynamic (RHD) simulation of supercritical accretion flows (Ohsuga et al. 2005, 
hereafter OMNM05), in which the viscous stress is expressed with the ap prescription, similar 
to our steady model. In Section 2, we describe the simulation model and steady model used in 
this paper. In Section 3, we check the self-similar assumptions by analysing the time-averaged 
simulation data in quasi-steady state. In Section 4, we compare the steady model solutions with 
the simulation results. Discussions are presented in Section 5 and we conclude with a summary in 
Section 6. 


2. Models 

In this section we briefly describe the simulation model and the steady model used in this 
paper. More details of the equations and numerical methods can be found in respective papers 
(OMNM05 & JWll). 


2.1. The simulation model 

The simulation is performed with the two-dimensional RHD code developed by Ohsuga et 
al.(2005). The calculation is carried out in spherical coordinates (r, 6 , (p), with a non-rotating black 
hole at the origin. The accretion flow is assumed to have axisymmetry (i.e. dld(f> = 0) and 
reflection symmetry relative to the equatorial plane {6 = n/l). Pseudo-Newtonian potential is 
adopted, which is given by i/f = -GM/(r - Vg) (Paczyhski & Wiita 1980), in which M is the black 
hole mass and Vg = IGMjc^ is the Schwarzschild radius. Radiative transfer is described by the 
flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981). It is assumed that 
the r^-component of the viscous tensor, tr^, is dominant. The dynamical viscosity coefficient is 
given by 

Pgas "I" flF'rad ,, 

V = a ---, ( 1 ) 

iZK 

where a is the viscosity parameter, Qk is the Keplerian angular speed, A is the flux limiter and 
^rad is the radiation energy density. It is basically the same as the a prescription of the viscosity 
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(Shakura & Sunyaev 1973), because A is almost 1/3 in the optically thick region, and the results 
show that for supercritical accretion flows the majority of the flow in the converged region is 
optically thick. 

The computational domain is set to be spherical shells of 3rg < r < SOOr^ and 0 < 6 < 
O.Sn, and is divided into 96 x 96 grid cells. The calculation is started with a hot, rarefied and 
optically thin atmosphere with no initial cold dense disk, and mass is injected continuously into the 
computational domain through the outer boundary (r = SOOr^, 0.457r < 9 < O.Sn) at a steady rate 
^input- The specific angular momentum of the injected mass is set to be the same as the Keplerian 
angular momentum at r = lOOr^. In this paper, two sets of simulation results are analyzed, both of 
which are taken from OMNM05 and have reached quasi-steady state for 3r^ < r < lOOr^. For both 
simulation runs, the parameters are set as M = IOMq, a = 0.1, y = 5/3, /i = 0.5 and Z = Z©, where 
y. Id and Z are the heat capacity ratio, mean molecular weight and metallicity, respectively. The 
mass injection rate is set to be 1000 Le/c^ and 3000 Le/c^ respectively, where Le is the Eddington 
luminosity. 


2.2. The steady model 

The steady accretion disk model with outflows is basically the same as presented in JW11. The 
calculation is carried out in spherical coordinates (r, 6, (p) with the gravitational center at the origin, 
and the flow is assumed to be steady(5/5t = 0) and axisymmetric(5/50 = 0). The Newtonian 
gravitational potential, O = -GM/r is adopted, and the viscosity is described by the a prescription 
tr^ = -ap, where p = Pgas + Pmd is the total pressure (magnetic pressure is not considered). It 
is assumed that the r^-component of the viscous tensor, is dominant. The energy equation is 
described by the advective factor, / = 2adv/2vis> so that a fraction / of the dissipated energy is 
advected as stored entropy and a fraction (1 - /) is lost due to radiation. Note that Newtonian 
potential, rather than the pseudo-Newtonian potential, is adopted here, so later in the comparison 
part, we will focus on the region not too close to the central black hole (r > lOr^), where relativistic 
effects can be neglected. 

Self-similar assumptions are adopted in the radial direction (Narayan & Yi 1995a; Xue & 
Wang 2005; etc.): 


pie)r-\ 

(2) 

Igm 


Vr(0) a/ 

(3) 

V r 


Igm 


Ve(9)^J ^ , 

(4) 
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p = p{e)GMr-’^-\ 


(5) 

( 6 ) 


and the hydrodynamic equations can be reduced to a set of ordinary differential equations (ODEs) 
about the variable 6 . This set of ODEs can be numerically solved with the symmetric boundary 
conditions relative to the equatorial plane if pinjl) is set to be 1, normalized by a scale factor if 
the effective accretion rate at a certain radius is set (Narayan & Yi 1995a; Xue & Wang 2005; etc.). 
Eour input parameters (Q',/,yequ>«) are required to calculate a solution, in which yequ is defined as 

y - 1 

yenu =-1" (7) 

where /3 = Pgas/p is the gas pressure ratio, yequ is defined to incorporate the influences of both gas 
and radiation pressure into the model. Eor pure ionized hydrogen (y = 5/3), in the case of extreme 
gas pressure domination (yS ^ 1), yequ = 5/3, while in the case of extreme radiation pressure 
domination (yS ^ 0), yequ = 4/3. 


The self-similar steady model here cannot describe the accretion flow structure in the whole 
space. Theoretically, this is because that, with the self-similar assumptions, the inflow accretion 
rate and outflow accretion rate scale with radius in the same way (both are proportional to ^“"). In 
the steady state, the total accretion rate should remain constant, which does not scale with radius. 
If the whole space can be described by the self-similar model, then there must be n=1.5 and both 
inflow and outflow accretion rate have to be constant. This means that there is no outflow (as 
outflow equals the difference between inflow accretion rates at different radii). In fact the n=1.5 
case has been solved by Narayan & Yi (1995a), in which all the streamlines are straight lines 
pointing at the central accretor and no outflow exists. So the self-similar steady model which 
contains outflow must have some boundary beyond which it cannot describe. Numerically, the 
calculation starts from the equatorial plane {6 = nil) and moves towards the polar axis {6 = 0). 
Both p and p decreases as 6 decreases, and at some polar angle they will get very close to 0 
simultaneously. If we continue the calculation beyond this angle, we will encounter numerical 
errors, so we take this polar angle 6 as the upper boundary of the steady model. 


It is worth mentioning that out of the four input parameters (a,/,yequ>^)> / does not appear in 
any differential terms in the equations, and thus can actually be variant in the 6 direction. The other 
three parameters are included in differential terms (see JW11 for detailed equations), but can also 
be set as functions of 6 , as long as the velocities, density and pressure do not change abruptly in 
space in the solution. To compare the steady model results with the simulation results, we obtain 
the input parameters (Qr,/,yequ,n) from the simulation data, which will be described in detail later. 


Here we adopt the self-similar assumptions presented in Narayan & Yi (1995a) and Xue & 
Wang (2005). It should be noted that in Narayan & Yi (1995a), the self-similar form of is 
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actually 

= v^( 0 ) sin e ( 8 ) 

although in their paper it is stated otherwise (possibly a typo; also see Tanaka & Menou 2006). If 
we use Eq. ([ 8 ]) to substitute Eq. (5), and leave other self-similar forms unehanged, we ean write 
another version of the eode. We run both versions of the eode and find that with the same input 
parameters (a,f,yequ,n), the results are aetually the same. It is not surprising beeause Eq.(5) is in 
faet equivalent to Eq.® if we inelude the sin 6 term into the v^( 6 ) term. 


3. Checking the self-similar assumptions 

In this seetion we eompare the self-similar assumptions used in the steady model with the 
simulation data. Erom Eqs.(2)-(6), we ean get the following relations (for brevity in this paper we 


use Tg’ to represent the logarithm with base 10): 

\gVr = ~\gr + ci{e), (9) 

\gve =—Xgr + C2{e), (10) 

Igv^ =—Xgr + c^ie), ( 11 ) 

\gp =-n\gr+ €^{6), (12) 

\gp =-{n+\)\gr + Ci{e). (13) 


Here Ci, C 2 , C 3 , C 4 and C 5 are dependant on the polar angle 6 , while they should remain eonstant for a 
fixed 6 under the self-similar assumptions for a eertain steady model solution. In JW11, we assume 
that n is eonstant, whieh means that for these relations the slopes should not ehange at different 
6 . The simulation data, on the other hand, do not adopt these assumptions, and do not neeessarily 
follow these relations. Here we fit the simulation data with linear models to eheek whether the 
self-similar assumptions in radial direetion are in good agreement with the simulation results. The 
data are averaged over t = 185 - 255 s of the simulation for Minput = 1000 Le/c^ and t = 182 - 252s 
for Minput = 3000Lb/c^, during whieh the simulations are in quasi-steady state, respeetively. As 
we are interested in the eonverged part of the simulation data exeept for the part very elose to the 
blaek hole, we foeus on the region from lOr^ to lOOrg in the simulation data. This eorrespond to 
the 23rd grid point to the 67th grid point. All the 96 rows divided in 6 direetion will be taken into 
eonsideration. 

The detailed fitting results and some diseussion are presented in the following subseetions. To 
evaluate the eurve fitting results, we make visual examination of the fitted eurves, as well as foeus 
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on two quantities: the fitted slopes and of fits. The fitted slopes are constant for velocities and 
directly connected with the input parameter n for density and pressure, and should remain constant 
for different 6 for a certain fit, as discussed above. is a goodness-of-fit statistic which is defined 
as 

ZCvi - 


R^ 


1 - 




(14) 


where y, is the value of a physical quantity from the simulation data, % is the fitted value from the 
linear regression, and y is the mean of y,-. R? measures how well the regression line approximates 
the simulation data points(the closer to 1, the better; an R? of 1 indicates that the regression line 
perfectly fits the data.). Table 1 summarizes all the fitting results. 


3.1. Fitting velocities 

The fitting results of v^, and are presented in this subsection. Linear fits have been 
performed at all inclination angles on the simulation grid, except for vq on the equatorial plane 
where ve = 0 due to reflection symmetry of the accretion flow. In this paper we define m = 
Miaputl(LElc^) for brevity. All the physical quantities are in cgs units. 

Figure 1 shows the linear fits of Ig Iv,-!, Ig Iv^l and Ig at polar angle 6 = 1.2732 rad in the 
inflow region, and Figure 2 shows the fits at 6 = n/A in the outflow region. We fit the absolute 
values of and Vg here, because in the simulation data there exist circulation patterns in the velocity 
fields in the rO plane, which is due to convection (see OMNM05). Figure 3 gives a general view of 
the fitting results of Ig |v^|, Ig |ve| and Ig at different polar angles. The upper panel corresponds 
to the R^ values, which are better for values closer to 1. The lower panel corresponds to the fitted 
slopes, which are better for values deviating less from each other on the same curve. For the fits 
of Ig Iv^l, we ignore 6 values for which there exists Vr = 0, and the left branch corresponds to the 
outflow region, while the right branch corresponds to the inflow region. 

The fits of Ig Iv^l displays two different types of behavior in the inflow and outflow regions, 
as shown in Figures 1 and 2. In the inflow region, the value of generally increases as radius 
decreases, but is strongly affected by the circulation patterns. Even after we remove the fits of 6 
values where Vr changes sign due to circulation patterns (i.e. the worst cases), the average R^ in 
the inflow region is only 0.31 and 0.59 for m = 1000 and m = 3000, respectively, with the best 
R^ valules of 0.64 and 0.91, respectively, as shown in Table 1. In the outflow region, the value of 
Vr first increases as radius decreases, then at some point it starts to decrease, as shown in the top 
panels of Figure 2, so that the fitted slopes of Ig Iv^l are mostly positive in the outflow region in 
Figure 3. This is because when the accretion flow gets close to the black hole, stronger gravity and 
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Table 1: Linear Fitting Results of Simulation Data (vs Ig r) 


Variables 




m = 1000 





Slope 

Mean Slope 

Ig Vr(inflow) 

0 .00- 

0.64 

0.31 

-0.94' 

-0.44 

-0.27 

Ig v^(outfiow) 

0.15 - 

0.80 

0.36 

0.07 - 

'0.88 

0.36 

Ig Vg 

0.21 - 

0.95 

0.81 

-1.34' 

-0.36 

-0.84 

Igv^ 

0.98 - 

0.99 

0.99 

-0.75 - 

' -0.56 

-0.65 

Igp 

0.87 - 

0.99 

0.94 

-1.43 - 

' -0.76 

-0.85 

Igp 

0.98 - 

0.99 

0.99 

-2.10- 

'-1.72 

-2.03 

Variables 




m = 3000 






Slope 

Mean Slope 

ig v^(infiow) 

0 .00- 

0.91 

0.59 

-1.20' 

-0.65 

-0.44 

ig v^(outfiow) 

0 .10- 

0.82 

0.31 

0.11 - 

'0.86 

0.38 

Ig Vg 

0.33 - 

0.98 

0.82 

-1.09' 

-0.24 

-0.74 

Igv^ 

0.97 - 

0.99 

0.99 

-0.73 - 

' -0.53 

-0.64 

Igp 

0.85 - 

0.99 

0.94 

-1.31 - 

' -0.55 

-0.74 

Igp 

0.98 - 

0.99 

0.99 

-2.02 - 

'-1.59 

-1.92 
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Fig. 1.— The linear fit eurves of Iglv^l, Ig |V(j| and Igv^ at polar angle 6 = 1.2732 rad in the 
inflow region. The black dots correspond to the simulation data, the solid lines correspond to the 
linear fitting, and the dashed lines correspond to the 95% confidence bounds. For m = 1000 and 
m = 3000, the of Ig \Vr\ fits are 0.6448 and 0.8857, respectively, and the fitted slopes of Ig \vr\ 
are -0.5981 (-0.7346, -0.4616) and -0.5775 (-0.6413, -0.5137), respectively, with 95% confidence 
bounds in brackets; the of Ig Iv^l fits are 0.7597 and 0.7979, respectively, and the fitted slopes of 
Igivel are -1.093 (-1.281, -0.9035) and -0.9937 (-1.147, -0.8399), respectively; the R^ of Igv^ fits 
are 0.9959 and 0.9952, respectively, and the fitted slopes of Igv^ are -0.6061 (-0.6181, -0.5941) 
and -0.5951 (-0.6078, -0.5824), respectively. 
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Fig. 2.— The linear fit eurves of Ig Iv^l, Ig |v6i| and Ig at polar angle 6 = n/A in the outflow region. 
The blaek dots eorrespond to the simulation data, the solid lines eorrespond to the linear fitting, and 
the dashed lines eorrespond to the 95% eonfidenee bounds. For m = 1000 and m = 3000, the of 
Ig \vr\ fits are 0.2478 and 0.3163, respeetively, and the fitted slopes of Ig \vr\ are 0.1399 (0.06495, 
0.2149) and 0.2064 (0.1131, 0.2998), respeetively, with 95% eonfidenee bounds in brackets; the 
of Igivel fits are 0.9209 and 0.982, respectively, and the fitted slopes of lg|ve| are -0.7764 (- 
0.8464, -0.7064) and -0.6071 (-0.6324, -0.5818), respectively; the R^ of Igv^ fits are 0.9995 and 
0.9999, respectively, and the fitted slopes of Ig are -0.698 (-0.7029, -0.6931) and -0.698 (-0.6999, 
-0.6961), respectively. 
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E? 



Slopes 



Fig. 3.— The fitted slopes and values of the linear fits of Ig Iv^l, Ig Ivel and Ig at different polar 
angles. The upper panel eorresponds to the values, and the lower panel eorresponds to the fitted 
slopes. Different markers represent different m and veloeity eomponents, as shown in the legend. 
For the fits of Ig Iv^l, we ignore 6 values for whieh there exists Vr = 0. 
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relativistic effects will gradually turn outflow into inflow (e.g. Figure 6 of OMNM05, in which 
inflow accretion rate becomes almost constant close to the black hole). Similar effects are also 
found in simulation works of ADAFs (Narayan et al. 2012; Yuan et al. 2012a), in which the inflow 
accretion rate are reported to be almost constant inside lOr^. This means that in the outflow region 
near the black hole, as radius decreases, Vr will decrease and eventually become negative, which 
corresponds to inflow. While this change usually happens inside lOr^, which we ignore in the 
linear regression, it does have influence on the profile of Vr in the radial direction, so that in the 
outflow region will start decreasing from some radius as r gets closer to lOr^. 

The fits of Ig are good, with fitted very close to 1 (0.99 on average for both m = 1000 
and m = 3000) and fitted slopes generally resembling the same value. So the self-similar model 
describes the radial distribution of well. The average fitted slopes of Ig are -0.65 and -0.64 
for m = 1000 and m = 3000, respectively, which differ from the slope of -0.5 in self-similar 
assumptions we use, but are still acceptable. The fits of Ig Iv^l are also good for most values of 6 , 
with an average of 0.81 and 0.82 for m = 1000 and m = 3000, respectively. The bad fits of 
Ig \v 0 \ appear near the equatorial plane and the polar axis, which is not surprising as in these places 
V 0 gets close to 0 due to symmetry and is thus much influenced by numerical errors. The fits of 
Ig Ivel are generally worse than those of Ig v^, as the profiles of are influenced by the circulation 
patterns we mentioned above, while those of are not. The slopes of the fits of Ig and Ig \vg\ are 
slightly different from -0.5, which means that their radial distributions are not strictly proportional 
to that of the corresponding Keplerian velocity. 

It should be noted that, near the equatorial plane close to the black hole, for both m = 1000 
and m = 3000, there is a small region which is severely influenced by the circulation patterns 
and displays another contour of Vr = 0, aside from the boundary between inflow and outflow 
region (c.f. the velocity field plot Figure 11 in Section 4). Although we ignore the 6 values for 
which there exists = 0, near this region the profiles are still much affected and the values 
of Vr drops significantly, making some fits of Ig Iv^l in the inflow region bad (with bad R^ values 
and positive slopes, as shown in Figure 3 and Table 1). If we regard the circulation patterns as 
perturbation caused by convection over a self-similar configuration, then it appears that the effect 
is much stronger in r direction than in 9 direction in spherical coordinates. 


3.2. Fitting density p 

The fitting results of density p are presented in this subsection. Figure 4 displays the linear 
fitting results of Igp at 0 = njl (i.e. on the equatorial plane) and 9 = 7r/4 for m = 1000 and 
m = 3000, respectively. All the physical quantities are in cgs units. The fitting results are good 
with R? values of 0.9591 and 0.9347 at 9 = njl, and 0.8902 and 0.9209 at 0 = njA, for m = 1000 
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and m = 3000, respectively. 



7.5 8 8.5 7.5 8 8.5 

lg(r) [cm] 


Fig. 4.— The fitting results of Ig p. The black dots correspond to the simulation data, the solid lines 
correspond to the linear fitting, and the dashed lines correspond to the 95% confidence bounds. For 
m = 1000, the of Igp fits are 0.9591 at 0 = njl, and 0.8902 at 0 = n/A, respectively, and the 
corresponding fitted slopes are -1.428 (-1.518, -1.338) and -0.7812 (-0.8656, -0.6968), respectively, 
with 95% confidence bounds in brackets. For m = 3000, the are 0.9347 at 0 = njl, and 0.9209 
at 9 = n/A, and the corresponding fitted slopes are -1.313 (-1.419, -1.206) and -0.7655 (-0.8345, 
-0.6965), respectively. 


Figure 5 gives a general view of the fitting results of density p at different polar angles. The 
upper panel corresponds to R^ values, which are close to 1 (0.94 on average for both m = 1000 
and m = 3000), so the fitting results are good. The lower panel corresponds to the fitted slopes at 
different polar angles, which generally resemble the same value on each curve, so the self-similar 
model describes the radial distribution of density in the the simulation data well. The average slope 
here is -0.85 for m = 1000, and -0.74 for m = 3000, which correspond to n ~ 0.85 and n ~ 0.74, 
respectively, according to Eq.(12). 
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Slopes 



Fig. 5.— The fitted slopes and values of the linear fits of Igp at different polar angles. The 
upper panel eorresponds to the values, and the lower panel eorresponds to the fitted slopes. 
Different markers eorrespond to different m, as shown in the legend. 
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3.3. Fitting pressure p 

The total pressure in the simulations eomes from both gas pressure and radiation pressure. 
The radiation stress tensor is almost isotropie in the optieally thiek region, so we have 


Pvad — Fjaij/S, 


(15) 


in whieh E^-ad is the radiation energy density (per unit volume). Note that for the simulation results, 
in the region very elose to the polar axis, the optieal depth eould be smaller than 1. The optieal 
depth T ean be approximated by "R (Kato et al. 2008): 


R = 


|V£| 

(^abs “1“ ^sca)pF 


1 

T 


(16) 


Figure 6 shows the eontour plots of R in the eonverged regions of the simulations for m = 1000 
and m = 3000. It ean be seen that most of the spaee is oeeupied by optieally thiek aeeretion flows 
for both simulations. The optieally thin regions are very small, and only have negligible influenee 
here. So the total pressure ean be written as 




Fig. 6.— Contour plots of R in the eonverged regions of the simulation data. R ~ 1/t indieates the 
loeal optieal depth. The left panel is for the simulation with m = 1000, while the right panel is for 
the simulation with m = 3000. 


P — Prad Pgas ~ ■^rad/3 + (y l)Fgas, (12) 

in whieh y is the heat eapaeity ratio whieh is taken as 5/3 in both simulations, and Eg^s is the gas 
internal energy density (per unit volume). 

Figure 7 displays the linear fitting results of Igp at 0 = njl (i.e. on the equatorial plane) and 
6 = njA for m = 1000 and m = 3000, respeetively. All the physieal quantities are in egs units. The 
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fitting results are good with E} values of 0.9997 and 0.9997 at 0 = nil, and 0.9865 and 0.9895 at 
6 = ;r/4, for m = 1000 and m = 3000 respeetively. 



7.5 8 8.5 7.5 8 8.5 

Ig(r-) [cm] 


Fig. 7.— The fitting results of Ig The blaek dots eorrespond to the simulation data, the solid lines 
correspond to the linear fitting, and the dashed lines correspond to the 95% confidence bounds. For 
m = 1000, the of Ig p fits are 0.9997 at 6 = n/l, and 0.9865 at 6 = n/A, respectively, and the 
corresponding fitted slopes are -2.105 (-2.116, -2.093) and -2.06 (-2.134, -1.985), respectively, 
with 95% confidence bounds in brackets. For m = 3000, the are 0.9997 atO = njl, and 0.9895 
at 6 = n/A, and the corresponding fitted slopes are -2.018 (-2.028, -2.007) and -1.901 (-1.961, 
-1.841), respectively. 


Figure 8 gives a general view of the fitting results of total pressure p at different polar angles. 
The upper panel corresponds to R^ values, which are very close to 1 (0.99 on average for both 
m = 1000 and m = 3000), so the fitting results are good. The lower panel corresponds to the 
fitted slopes at different polar angles, which generally resemble the same value on each curve, so 
the self-similar model describes the radial distribution of total pressure in the simulation data well. 
The average slope here is -2.03 for m = 1000, and -1.91 for m = 3000, which corresponds to 
n ~ 1.03 and n ~ 0.91, respectively, according to Eq.(13). Note that the values of n obtained here 
do not agree completely with the n values obtained from the fits of p, but are quite close. This will 
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be discussed in detail later in the paper. 


4. Comparing the steady model solutions with simulation results 

In this section we compare the steady accretion disk model solutions with the simulation 
results. First we determine the input parameters for the steady model from the simulation data. 
Then the steady model solutions corresponding to these parameters are calculated and compared 
with the simulation results. 


4.1. Determining input parameters 


Four input parameters are required for the steady accretion disk model, namely (a,f,yequ,n). 
The simulation data are calculated with or = 0.1 and y = 5/3 as mentioned in section 2. According 
to Eq.(7), 7equ is determined by both y and /3. Figure 9 displays the contour maps of gas pressure 
ratio jS obtained from the simulation data. It can be seen that, in the converged region of the simula¬ 
tions, jS is generally below 0.2, except for a very small region close to lOr^. The converged regions 
in both simulations are radiation-pressure dominated, with a maximum/J^ax = 0.386 (correspond¬ 
ing to yequ = 1.41) for m = 1000, and = 0.344 (corresponding to yequ = 1-40) for m = 3000. 
So here we can safely take yequ = 4/3, which corresponds to extremely radiation-pressure domi¬ 
nated accretion flows. 


The parameter n can be obtained from either the fits of density p or the fits of total pressure p, 
as indicated by Eqs.(12) and (13). As discussed above, we get n ~ 0.85 and n ~ 1.03, from the fits 
of p and p, respectively, for m = 1000. Eor m = 3000, we get n ~ 0.74 and n ~ 0.91, respectively. 
However, the relation indicated by Eq.(13) is based on the assumption: 


oc 


GM 


(18) 


r 

in which is the sound speed. As discussed in section 2, the profile of actually scales in 
the radial direction with on average for m = 1000, and on average for m = 3000. 
It is possible that Cs also deviates from the scaling in the radial direction. So it is better to 
determine n with the density profiles. Here we set n = 0.85 for m = 1000, and n = 0.74 for 
m = 3000, from the fits of density p. It does not differ a lot from the values obtained from the fits 
of pressure p, anyway. 


The energy equations adopted in the simulations are 


BE 

gas 


-|- V • (Egas^^) Pgas^ ■ ^ — Q pj “1“ PC/CabsEradj 


dt 


(19) 
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Fig. 8.— The fitted slopes and values of the linear fits of Ig p at different polar angles. The 
solid curves correspond to the values, and the dashed curves correspond to the fitted slopes. 
Different markers correspond to different m, as shown in the legend. 
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Fig. 9.— Contour plots of gas pressure ratio jS in the converged regions of the simulation data. 
The left panel is for the simulation with m = 1000, while the right panel is for the simulation with 
m = 3000. 


— — + V • (£radV) + Vv : P = -V • F + P7 - pC/CabsF^rad, (20) 

ot 

in which v is the velocity, P is the radiation pressure tensor, q'^ is the viscous heating rate, j is 
the emissivity per unit mass, and /Cabs is the absorption opacity. Eq.(19) is the energy equation 
of the gas, and Eq.(20) is the energy equation of the radiation. The left hand side of these two 
equations, excluding the time derivative terms, are the advection of the gas internal energy and the 
advection of the radiation field, respectively. In the steady model, the advection term (^adv includes 
both sources of advection. Correspondingly, if we sum up Eq.(19) and (20), we can get 


diEgas Efad) 
dt 


^adv — Q V * F . 


( 21 ) 


V • F is similar to the radiation cooling in one-dimensional accretion disk models. It represents the 
change of the radiation energy in a fixed region due to the transportation via the radiative flux F. 
The advective factor / can be calculated from Eq.(21): 




^adv 


1 - 


^(£gas+grad) 

dt 


+ V-F 


( 22 ) 


As we focus on the converged region in the simulation which has achieved quasi-steady state, the 
time derivative terms are typically much smaller than other terms (10^^ - 10“^ times q'^), so the 
values of / mainly depend on q"" and V • F. We then average / over time and radial direction, 
to get its profile in 6 direction. Eigure 10 displays the result. Basically / is close to 1 on the 
equatorial plane (where 6 = njl), decreases as inclination decreases, and then increases again as 6 
gets close to the polar axis, and values of / can get very large (greater than 10, not shown in Eigure 
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10 to display more details for lower values of /) near the polar axis. This can be explained by the 
photon trapping effect (Begelman 1978; Takeuchi & Mineshige 2009; also see OMNM05) and the 
existence of strong outflow near the polar axis. The photons generated near the equatorial plane 
are more effectively trapped than those generated at higher latitude; as 9 decreases, it is easier for 
the generated photons to escape the accretion flow, so / decreases accordingly. However, below a 
certain value of 9, more photons are generated in the accretion flow at lower latitude and carried 
inside by the radiative flux than those escaping to higher latitude, so that V • F becomes negative, 
which causes / to become larger than 1. These photons are then carried outward by the strong 
outflow near the polar axis, together with kinetic and internal energy of the gas, in the form of <7adv 
This is also verified later as we investigate the kinetic energy flux carried by outflows in Section 
4.4. It should be noted that the region near the polar axis is not described in the steady model, so 
it does not influence our calculation. The dashed lines in Figure 10 represent the boundary of 9 in 
the steady model calculations. 

The values of / are comparably larger in the simulation for m = 3000 than those for m = 1000. 
It is natural because for supercritical accretion flows, larger accretion rate corresponds to heavier 
photon trapping, which causes / to be larger. 


4.2. Comparing the velocity fields 

In this subsection, we present the velocity field plots of the simulation results and the steady 
model solutions, as shown in Figure 11. The solid curves in each panel correspond to the surface 
of the inflow region at which Vr = 0. The dashed lines correspond to the upper boundary of the 
steady model calculations beyond which our self-similar model can no longer describe (see JW11 
for more details). In each velocity field plot of the simulation data, there is a small region wrapped 
by a solid curve near the equatorial plane. That corresponds to a circular pattern of the accretion 
flow, which arises from convection, as mentioned in Section 3. 

It can be seen that for both accretion rates, the simulation result displays a smaller region of 
inflow than the steady model result. For m = 1000, the upper surface of the inflow region resides 
at 0 = 1.17 rad (~ 67.0°) on average for the simulation, while it is at 0 = 0.972 rad (~ 55.7°) 
for the steady model. For m = 3000, the upper surface of the inflow region resides at 9 = 1.06 
rad (~ 60.7°) on average for the simulation, while it is at 0 = 0.986 rad (~ 56.5°) for the steady 
model. Besides that, in the outflow region of the simulation, the outward radial motion dominates 
the accretion flow, while in the steady model the flow follows a combination of and v^. This 
indicates that either Vr is much larger, or vg is much smaller, in the outflow region of the simulation 
results than those of the steady model results (in Section 4.3 we will see that Vr is underestimated 
at high latitude in the steady model). 
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9 [rad] 


Fig. 10.— Profiles of the advective factor f in 6 direction, calculated by taking its average over 
time and radial direction. Black dots correspond to the simulation data, the solid curves are data 
smoothed by a locally weighted linear least-squares regression method, and the dashed lines rep¬ 
resent the boundary of the steady model calculations in 6 direction. 
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Fig. 11.— The velocity field plots of the simulation (upper panels) and steady model (lower 
panels) results. In each panel, the solid curves correspond to the surface of the inflow region, while 
the dashed line corresponds to the calculation limit of the steady models. The solid curves near 
the equatorial plane in the upper panels correspond to circular patterns of accretion flow in the 
simulation results, which arise from convection. 
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As the velocity field plots of the simulations here are based on time-averaged data, the convec¬ 
tive motions in large circular pattern are not so obviously shown as in a snapshot of the simulation 
(cf. OMNM05). However, convection still works, though not dominant, as a transport mechanism 
of energy, mass and angular momentum. This will be discussed in Section 5. 


4.3. Comparing the variable profiles in 6 direction 


In this subsection, we compare the profiles of velocities, density and pressure along the 6 
direction, obtained from the steady model solutions and the simulation data. In the steady model, 
the profiles are presented in the form of Vr(6), Vg(6), v^{6),p{6) and p{6). The corresponding forms 
in the simulation results are 


^r,sim(^) — 

Psim(0) — 
Psim(^) — 


Vr 


Vk 

sim 

Vk 

V(/>,sim 

Vk 

Psim 

Psim,6=;r/2 
Psim 
PsimVj^ 


Psim(^) 


Psir 


1 


(23) 

(24) 

(25) 

(26) 

(27) 


Psim,6=;r/2 Vj^ 

in which the subscript ’sim’ indicates the data from simulation results and Vk is the Keplerian 
velocity 


Vk = 


GM 


r - r„ 


(28) 


All these calculations are done at the same radius. 


Figure 12 displays the profiles of v^, vg and along the 6 direction. Note that as discussed 
in Section 3, the simulation data do not strictly follow the self-similar assumptions, especially for 
Vr and Vg, so the profiles for the simulation results are also dependant on the radius. Here we 
present the profiles at three different radii lOr^, 49rg and 99rg (all these are on grid points of the 
simulation), distinguished by different markers. The steady model results are represented by the 
solid curves, which are not so smooth as presented in JW11, due to the fact that the calculations are 
carried out with variable values of / (see Section 4.1). The simulation model is based on pseudo- 
Newtonian gravity while the steady model is Newtonian, so Vk is slightly different from each other. 
However, as we focus on the region between 10 and 100 r^, the difference is negligible. 

As discussed in Section 3, (and partly v^) does not follow the self-similar assumptions in 
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Fig. 12.— The distribution of v^, Vg and in 0 direction, all normalized by the local Keplerian 
velocity. The profiles corresponding to the simulation data are displayed for three different radii of 
lOrg, A9rg and 99rg, indicated by black dots, green squares and cyan triangles, respectively (see the 
online version for the colored figure). The steady model results are represented by solid curves. 




























25- 


the simulation results, likely due to the eireular patterns arising from conveetion. However, the 
steady model profiles of and generally follow the same eurve shapes of those obtained from 
the simulation data, exeept for the profile of near its upper boundary, whieh arises from the fluc¬ 
tuations of the steady model calculations near the upper boundary, tends to be underestimated 
in the steady model for large radii, which corresponds to the difference in the flow motion of the 
velocity field plots mentioned in Section 4.2. 

From the fits of in Section 3, we know that follows the relation oc on average for 
m = 1000 and oc for m = 3000. That deviates from the radial scaling of Vk oo so that 
although the fitting results of have an average E} value of 0.99 for both m = 1000 and m = 3000, 
the profiles of at different radii do not coincide. The profiles of in the steady model results 
tend to increase as the polar angle decreases, which seems to not agree with the simulation results. 
However, the profile of angular velocity Q actually agrees with the simulation results, which will 
be shown in Section 5. 

Figure 13 displays the profiles of p and p in the 6 direction. The profiles obtained from the 
simulation data are shown at three different radii, for the same reason as discussed above. Note that 
although the radial distribution of total pressure p follows the self-similar assumptions well with 
an average R? value of 0.99 for both m = 1000 and m = 3000, as discussed in Section 3, the input 
parameter n for the steady model is chosen based on the fits of p, and the exponent is not the same 
(n ~ 0.85 and 0.74 for p while n ~ 1.03 and 0.91 for p, for m = 1000 and 3000, respectively). That 
is why the three profiles of p at different radii of the simulation do not coincide. It can be seen that 
the profiles from the steady model results generally follow the same curve shape as the simulation 
profiles. The pressure p near the equatorial plane are overestimated in the steady model results for 
both accretion rates. 

Generally speaking, the steady model results agree with the simulation data not only qual¬ 
itatively, but also quantitatively, deviating from the simulation data only several times at most. 
Considering that the steady model is based on very simple assumptions and calculations, the result 
is quite satisfying. 


4.4. Mass, momentum and energy fluxes in outflows 

In this subsection, we investigate the mass, momentum and kinetic energy fluxes in the out¬ 
flows, calculated from the steady model solutions and the simulation data, as shown in Figure 14. 
The solid and dashed curves correspond to the fluxes obtained from the simulation data and the 
steady model solutions, respectively, and the dotted lines correspond to the calculation boundaries 
in the steady model solutions. In the steady model, the density on the equatorial plane is set to be 
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6 [rad] 


Fig. 13.— The distribution of p (normalized to p(nl2)) and p (normalized by Eq.(23)) in 6 di- 
reetion. The profiles eorresponding to the simulation data are displayed for three different radii of 
lOr^, 49rg and 99rg, indieated by blaek dots, green squares and eyan triangles, respeetively (see the 
online version for the eolored figure). The steady model results are represented by solid eurves. 
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1, which can be regarded as normalized by a scale factor. To translate the fluxes obtained from the 
steady model into real physical units, we need to calculate this scale factor, which is obtained by 
setting the mass inflow rate at 500 (which is the outer boundary in the simulations) to be the same 
as the mass supply rate parameter in the corresponding simulation calculations. The momentum 
flux presented here considers the total momentum, including the ^ component. 



Fig. 14.— Mass, momentum and kinetic energy fluxes in the outflows. The solid and dashed 
curves correspond to the fluxes obtained from the simulation data and the steady model solutions, 
respectively. Dotted lines indicate the calculation boundaries in the steady model solutions. 


The mass flux obtained from the simulation data peaks at 0 = 0.8914 rad (~ 51.1°) for 
m = 1000, and 6 = 0.5381 rad (~ 30.8°) for m = 3000, respectively. Mass flux also remains 
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high around respective peaks, which means that the outflowing mass is overall directed toward the 
polar angle around the peak. At the outer boundary of the converged region (r = lOOr^), ~94% of 
the mass inflow are driven away as outflow. The mass flux obtained from the steady model peaks 
at 0 = 0.8382 rad (~ 48.0°) for m = 1000 and 6 = 0.7852 rad (~ 43.4°) for m = 3000, respectively. 
Both peaks fall in the main outflow region as shown in Fig. 13. There are still significant outflows 
beyond the calculation boundary of the steady model, and the mass outflow in the calculated region 
of the steady model only takes up around 9% of the mass inflow. 

For the momentum and kinetic energy fluxes, although those obtained from the simulation 
data basically agree with those obtained from the steady model in the regions that have been 
calculated, they keep increasing until near the polar axis, reaching maximum values on the 2nd 
calculation grid point (9 = 0.1749 rad, ~ 10.0°) for mass flux, and on the 1st calculation grid point 
(6 = 0.0724 rad, ~ 4.1°) for kinetic energy flux. Significant momentum and kinetic energy fluxes 
exist in the regions around the polar axis that are not calculated by the steady model, and should 
be taken into consideration if one applies our steady model to observations. 

As the mass supply rate increases from m = 1000 to m = 3000, the fluxes appear to increase 
more for the steady model solutions than for the simulation results. This is likely because that the 
simulation results have not achieved converged state from lOOr^ to SOOr^, so that increasing the 
mass supply rate at the outer boundary has less impact on the outflows at lOOr^ in the simulation 
results. 


5. Discussion 

5.1. Validity of self-similarity and steady solutions 

The steady model we adopt here is based on the assumption of self-similarity. Physically, 
the assumption is based on the idea that in the problem, the only length scale of interest is r, and 
the only frequency is Qk, so that components of velocity, as well as the sound speed, should scale 
with radius as rflx (Narayan & Yi 1995a). The radial distribution of p is not quite clear, so it is 
described with a variable n in the form of p oc r“”. Due to the relation that the isothermal sound 
speed Cs = yjplp, we can get the radial distribution of p as p cc In the extreme case that vg is 
assumed to be 0, the value of n can be calculated to be 1.5 via the continuity equation. In this case, 
all the stream lines of the accreted material would be straight lines pointing directly at the central 
accretor, and the accretion flow would be composed of pure inflow. For n < 1.5, outflows will be 
generated from the accretion disk, n > 1.5 would indicate that the disk is supplied with mass in the 
6 direction as ’inflow wind’, which is unlikely to happen in real cases. Self-similar assumptions 
have been widely adopted in steady accretion disk models (e.g. Narayan & Yi 1994, 1995a; BB99; 
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BB04; Xue & Wang 2005; Gu et al. 2009; JWll; Begelman 2012; Gu 2012, 2015; etc.). However, 
the validity of these assumptions are not quite clear. 

In Section 3, we have checked these assumptions with two samples of RHD simulations of su¬ 
percritical accretion flows. The radial distribution of is close to the self-similar form, especially 
near the equatorial plane, with the exponent of r around -0.6 (Figure 3). The radial distribution 
of Vr and ve, however, do not follow the self-similar forms so well. This is mostly due to the 
circular patterns in the velocity fields of the simulations. The self-similar forms require that 
and remain in the same direction for the same polar angle 6, therefore the stream lines can not 
form enclosed rings in the self-similar region, which is obviously violated by the circular patterns. 
These patterns are not obvious in the time-averaged velocity field plots, but can be seen clearly in 
a snapshot figure (such as Figure 4 in OMNM05). The radial distribution of in the outflows is 
also influenced by the strong gravitational force near the central black hole, which causes the value 
of Vr in the outflows to drop significantly near lOr^. Density p and pressure p are in agreement 
with the relations p oc and p cc as we assume. In our steady model, the mass inflow 
rate Mintiow ^ The parameter n obtained from the simulation data is smaller for m = 3000 

than m = 1000, which implies that stronger outflows are generated when the mass supply rate is 
increased. 

The distribution of velocities, density and pressure in 6 direction calculated by our steady 
model is basically in agreement with the simulation result, displaying only serval times of differ¬ 
ence quantitatively at maximum (see Section 4.3). The largest deviations appear for the 9 profiles 
of azimuthal velocity and pressure p. For v^, the curve shape seems different from the simu¬ 
lation results in Figure 12. Actually it is not that different, if we show the 6 profiles of angular 
frequency G instead (Figure 15). In Figure 15, G/Gk decreases as 6 decreases for both the steady 
model solutions and the simulation results, which is natural as ’isorotes’ (surfaces of constant G) 
tend to be less curved than spheres. Note that although both vk and Gk are set as constants for the 
same spherical radius, the rotational motion is actually relative to the cylindrical radius, so that the 
curve shape of Q/Qk is different from that of v<^/vk. Pressure p is overestimated at low latitude and 
underestimated at high latitude in our model, compared with the simulation data. We conjecture 
that this arises from the fact that convective motions in the simulation act as an additional mecha¬ 
nism of transporting energy outward, reducing the energy density and consequently the pressure at 
low latitude and increasing them at high latitude. The convective motion may also play some role 
in transporting angular momentum, but this is not as clear as the pressure distribution. 

It should be noted that the analyses in this paper are based on two samples of RHD simulation 
runs, and need further confirmation from more simulation data. 2D radiation-magnetohydrodynamic 
(RMHD) simulations (Ohsuga et al. 2009; Ohsuga & Mineshige 2011) and 2.5D special relativis¬ 
tic RMHD simulations (Takahashi & Ohsuga 2015) of supercritical accretion flows have been 
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Fig. 15.— The distribution of Q (normalized by eorresponding loeal in 0 direetion. The 
profiles eorresponding to the simulation data are displayed for three different radii of lOr^, 49rg 
and 99rg, indieated by blaek dots, green squares and eyan triangles, respeetively (see the online 
version for the eolored figure). The steady model results are represented by solid eurves. 

performed in reeent years, but these simulations are mueh more time-eonsuming and thus only 
have limited ealeulation time. The general inflow-outflow strueture in these simulations is similar 
to that of our steady model, with inflow near the equatorial plane and outflow above the inflow 
region. They also show new features sueh as magnetieally eollimated jets around the polar axis, 
whieh are not deseribed by the steady model. However, limited ealeulation time, restrieted eompu- 
tational domain and the influenee of the initial eonditions make the eurrent simulations not suitable 
for detailed comparison, and we leave this topic for future work. 


5.2. Convection 

In our steady model, it is assumed that the energy transport mechanism is mainly based on 
advection and radiation, while heat conduction and convection are neglected. While our model is 
self-consistent, there are certainly other possible solutions with alternative prescriptions. As dis¬ 
cussed in the above section, convection seems to have some impact on the structure of supercritical 
accretion flows. The convective transport is essentially non-local, which acts in large circulation 
patterns in the snapshots of the simulation data. However, the time-averaged data do not display 
obvious circulation patterns (see Figure 11). Even so, the effects of momentum and energy trans- 
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port through convective motions still exist, and need to be considered to improve the agreement 
between steady solutions and the simulation results, which to some extent represent real astrophys- 
ical accretion flows. 

In literature, there is a series of self-similar models for accretion flows with convection: the 
’adiabatic inflow-outflow solutions’ (ADIOS; BB99, BB04; revised in Begelman 2012, which no 
longer requires convection to be the dominant mechanism of energy transport). The BB99 version 
of ADIOS is a one-dimensional (ID), height-integrated, radially self-similar model of steady ra- 
diatively inefficient accretion flows. It is basically a variation of the ID self-similar ADAF model 
proposed by Narayan & Yi (1994), with a variable accretion rate m oc {p correspond to 1.5 - n 
in our model). This model is expanded and redefined in BB04. They propose that the model is 
applicable to both ’overfed’ (i.e. supercritical) and ’underfed’ scenarios of accretion disks and the 
model is expanded to 2D. The 2D ADIOS model also contains outflow, and poloidal flow in the 
model is quadrupolar, inward at low latitude and outward at high latitude, which is similar to our 
model. In their model, convection is included as the main mechanism of energy transport. How¬ 
ever, this is not reflected in the energy equation, which is described with polytropic relation (they 
only require p oc'^ when an element of gas changes density) and the energy conservation equation 
is used to solve for the convective energy flux after the structure of the flow has been obtained. 
The main postulation about convection in BB04 is that the convective motion is performed along 
the gyrentropes (surfaces of constant specific angular momentum, Bernoulli function and entropy, 
which coincide with each other when the disc is marginally stable to the second Hpiland criterion), 
which does not agree with the simulation results here, as convective motions in the simulation data 
actually operate in large circular patterns. 

While both our model and the ADIOS model are based on self-consistent assumptions, our 
model is more suitable to be applied to the analysis of the simulation data from OMNM05 mainly 
for two reasons. The first is that, in ADIOS, the energy transport due to radiative flux is assumed 
to be zero, which is not applicable to the simulations of supercritical accretion flows here (cf. 
Section 4.1). The transport of energy is mainly due to advection and radiation, which agrees with 
the assumptions of our steady model. More importantly, in both our model and the simulation 
data, viscosity still works in the outflow region, while in the ADIOS model it is assumed to no 
longer work in the outflow part. The outflow-disk model described in BB04 calculates the disk 
part and outflow part separately. The outflow is assumed to launch from a ’thermal front’ where 
the convective motions quickly dissipate, increasing the entropy of the gas. It is also assumed that 
in the outflow part, the viscous torque no longer act and the viscous transport of angular momentum 
and energy also stops. This is quite different from our modef, as we calculate the structure of the 
whole accretion flow together, and the outflow is identified after the structure has been obtained. 
Obviousiy our model agrees more with the simulations in this sense. We are also planning to 
include convection in our steady model in future. 
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The ADIOS model is revised by Begelman (2012), in whieh the outflow ealeulation is mueh 
improved, as the outflow strueture is no longer assumed to be laminar and inviseid. However, the 
vertieal strueture of aeeretion flow is not obtained in the revised ADIOS model, so it is diffleult to 
eompare with our results here. It is suggested in the revised ADIOS model that the mass flux index 
for supereritieal aeeretion should be less than 1, which corresponds to n > 0.5 in our model (note 
that we use different parameter notifications). This agrees with our results in this paper. 


5.3. The ’underfed’ case 

As mentioned in introduction, outflows are likely to be generated in both the ’overfed’ and the 
’underfed’ case. We have discussed the ’overfed’ case with our steady model and simulation data. 
However, our steady model is established as a general model which can describe both the ’overfed’ 
and the ’underfed’ case by adjusting the input parameters. It is also an interesting topic to compare 
our steady model with the simulation data in the ’underfed’ case, which we are planning to do in 
future work. Here we just give some general notes in the ’underfed’ case. 

The ’underfed’ case corresponds to a hot, optically thin accretion flow which is radiatively 
inefficient and thus has an advection-dominated energy transport. Therefore the advective factor 
can be taken as 1. Radiation pressure is also very small so that ygqu can be taken as 5/3, which 
corresponds to gas-pressure dominated accretion (note that some simulations take into account the 
magnetic pressure, in which case yequ should be adjusted accordingly). This leaves a and n to 
be obtained from simulation data, while the validity of self-similar assumptions should also be 
checked. 

The viscous parameter a can either be obtained from the form of shear stress adopted in 
hydrodynamical simulations (e.g. Yuan et al. 2012a), or be calculated from simulation data of 
MHD simulations (e.g. Narayan et al. 2012). The parameter n can be calculated from radial 
density profile of the simulation data. For example. Yuan et al. (2012a) obtained p oc 
p oc oc oc in the case of a = 0.001, and p oc p oc with velocity 

index around -0.5 (observed from their Figure 4) in the case of a = 0.01, both for r > lOr^. 
Narayan et al. (2012) also find oc for r > lOr^. It appears that there are positive evidence 
for the validity of self-similar assumptions in the ’underfed’ case, although detailed investigation 
still awaits to be made. There are also debates over the importance of convection in the ’underfed’ 
case (e.g. Yuan & Bu 2010; Narayan et al. 2012; Yuan et al. 2012a, 2012b), which will be 
investigated in future work. 



6. Summary 


We make comparison between our steady accretion disk model containing outflows and two 
samples of 2D RHD simulation of supercritical accretion flows. The steady model is based on 
radial self-similar assumptions of velocity, density and pressure, which are checked with the sim¬ 
ulation data. In the converged region of the simulation data, excluding the part too close to the 
central black hole, azimuthal velocity v^, density p and total pressure p basically follow the self¬ 
similar assumptions, while radial velocity Vr does not, which is likely due to the circular pattern of 
accretion flow in the rO plane of the simulation, caused by convection. The radial distribution of Vr 
in the outflow region is also influenced by the strong gravitational force near the central black hole, 
which causes the value of in the outflows to drop significantly near lOr^. Polar velocity Ve some¬ 
what follows the self-similar assumptions, although it is also influenced by convection. Physically, 
convection acts as an additional mechanism of transporting momentum other than those consid¬ 
ered in our self-similar model, so that profiles of Vr and ve are disturbed and deviate from their 
self-similar forms. The fact that Ve follows self-similarity better than implies that the effects of 
convection are stronger in the radial direction, giving us some hint on how to treat convection in 
future work. 

Then we calculate the solutions of the steady model, with input parameters based on the sim¬ 
ulation data. In the steady solutions, the distribution of physical quantities in 6 direction basically 
agree with the simulation results. The agreement is good not only qualitatively, but also quantita¬ 
tively, as the steady model results deviate only several times from the simulation results at most. 
The result of comparison is satisfying, considering that the steady model is based on very simple 
assumptions and calculations. In the simulation results, outflowing mass is overall directed toward 
polar angle of 0.8382 rad (~ 48.0°) for m = 1000, and 0.7852 rad (~ 43.4°) for m = 3000, and 
~94% of the mass inflow are driven away as outflow (at lOOr^), while outward momentum and 
kinetic energy fluxes are focused around the polar axis. There exist significant mass, momentum 
and kinetic energy fluxes in the regions that are not calculated by our steady model, and attention 
should be paid when the model is applied to other theoretical or observational studies. The radial 
velocity Vy is underestimated at high latitude, and the total pressure p is overestimated at low lati¬ 
tude in the steady model. We conjecture that if convection is included as an additional mechanism 
of transporting mass, energy and angular momentum, this disagreement may be alleviated. In the 
two samples of simulation data in this paper, convection is less important than advection and radi¬ 
ation, so our steady model still holds in principle. Convection may play a larger role in optically 
thin accretion flows (e.g. Yuan & Bu 2010, but see Narayan et al. 2012 for a different opinion), 
which will be investigated in future work. We are also planning to include convection in the steady 
model in future. 

The analyses in this paper are based on two samples of simulation runs, and need further 
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confirmation from more simulation data. As the first step, we compare our steady model with 
RHD simulation of supercritical accretion flows based on ap prescription of viscosity. As our 
steady model is parameterized and can be adjusted to correspond to both ’overfed’ and ’underfed’ 
accretion flows, we would like to apply our model to simulations of ’underfed’ accretion flows in 
future work. Besides that, there have been advanced MHD simulations recently in which viscosity 
is generated through magneto-rotational instability (MRI) self-consistently (e.g. Ohsuga et al. 
2009; Ohsuga & Mineshige 2011; Narayan et al. 2012; Yuan et al. 2012b; Takeuchi et al. 2013; 
etc.). It is not clear whether the results of these complicated simulations can be reproduced by 
simple analytic models, which will be a topic in our future work. 


We thank the referee for making very helpful comments and suggestions. 
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